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Pulsed ultrasonic techniques can be and 
have been used to examine the interface 
conditions of a bonded structure. To provide 
a theoretical basis for such testing tech- 
niques we model the structure as a layer on 
top of a half-space, both of different 
elastic properties, with various interface 
bonding conditions. The exact dynamic 
Green's tensor for such a structure is ex- 
plicitly derived from the three-dimen- 
sional equations of motion. The final solu- 
tion is a series. Each term of the series 
corresponds to a successive arrival of a 
"generalized ray" and each is a definite 
line integral along a fixed path which can 
be easily computed numerically. Willis' 
method is used in the derivation. A new 
scheme of automatic generation of the 
arrivals and ray paths using combinatorial 
analysis, along with the summation of the 
corresponding products of reflection coeffi- 
cients is presented. A FORTRAN code is 
developed for computation of the Green's 
tensor when both the source and the de- 
tector are located on the top surface. The 



Green's tensor is then used to simulate 
displacements due to pulsed ultrasonic 
point sources of known time waveform. 
Results show that the interface bonding 
conditions have a great influence on the 
transient displacements. For example, when 
the interface bonding conditions vary, 
some of the first few head waves and regu- 
lar reflected rays change polarities and 
amplitudes. This phenomenon can be used 
to infer the quality of the interface bond 
of materials in ultrasonic nondestructive 
evaluation. In addition the results are 
useful in the study of acoustic microscopy 
probes, coatings, and geo-exploration. 
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Introduction 



The dynamic Green's tensor of a structure is the fun- 
damental solution to the transient mechanical wave 
problem of the structure. The theoretical prediction of 
the behavior of transient waves for arbitrary extended 
source distributions in space and time can be obtained 
once the time-domain dynamic Green's tensor is known. 
The ability to predict the behavior of transient waves in 
a structure is important in the development and under- 
standing of nondestructive evaluation techniques using 

1 Deceased; was on leave from Institute of Acoustics, Academia 
Sinica, Beijing, People's Republic of China. 



ultrasonics or acoustic emission and in other problems 
of wave propagation in solid and liquid media. 

In recent years, the "generalized ray" expansion tech- 
nique has been applied to compute the Green's tensor in 
a solid half-space or a solid infinite plate. The governing 
differential equation of motion is transferred to the 
Fourier or Laplace domain and the solution in the form 
of a series is obtained algebraically. Basically, there are 
two methods to invert the series from the transformed 
domain to the time-space domain. One is the well 
known Cagniard-de Hoop inversion method [1]; the 
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other is a method developed by Willis in 1973 [2], 
Willis' method uses the Fourier transform and expresses 
the resulting transform as a series of "generalized rays" 
and then inverts the series term by term. For three-di- 
mensional problems, as in the Cagniard-de Hoop 
method, only one integration remains. The integration 
path is always around a unit circle and is therefore 
"fixed" to some extent, but explicit evaluation of the 
integrand requires the numerical solution of an algebraic 
equation for each integration variable. Application of the 
Cagniard-de Hoop method requires a detailed discus- 
sion of the structure of a moderately complicated alge- 
braic function accompanying the transform of the inte- 
gration path. However, the integration path of the Willis 
inversion method is "fixed" and succeeds in avoiding the 
explicit discussion of the structure of the algebraic func- 
tion and so is applied rather easily even in the an- 
isotropic case. This is the main advantage of the Willis 
inversion method. The basis for carrying out the Willis 
inversion is that the solutions of elastodynamic problems 
are homogeneous functions of time, t, and position, x; 
that is, the problems are self-similar [2]. In addition, the 
displacements rather than the potentials are used, so the 
derivations are simplified, and the derivatives of the 
Green's tensor about spatial coordinates are easily ob- 
tained as well. These derivatives can be interpreted as 
the displacements due to dipole sources. The three-di- 
mensional transient Green's tensor for an isotropic plate 
and the two-dimensional transient Green's tensor for an 
anisotropic layer on an isotropic half- space were solved 
using Willis' inversion method [3,4]. 

In this paper formulae are derived for computing the 
three-dimensional transient Green's tensor for an 
isotropic layer overlay on an isotropic half-space. To 
understand the influence of the interface bonding condi- 
tion on the behavior of transient waves, a welded inter- 
face, a liquid coupled interface, and a "vacuum" inter- 
face are considered. The results for the case of the liquid 
coupled interface are obtained by artificially casting the 
boundary conditions into a matrix form similar to that 
for the case of the welded interface. The results for the 
"vacuum" interface are obtained by considering the 
layer with no half-space. The last case has been com- 
puted and experimentally confirmed previously, thus it 
can be checked with independent results. 

There are many rays which arrive at the observation 
point (detector) at the same time owing to the multiple 
reflection and the mode conversion of the incident P 
(longitudinal) ray or S (shear) ray emitted from the force 
source in the layer. These rays are kinematically equiva- 
lent and are called "kinematic analogs". Obviously, it is 
not necessary to separately compute the contribution of 
each ray to the integration. The question of how many 
kinematically equivalent rays arrive at the detector at the 



same time for a given configuration is a problem of 
combinatorics. It is quite a complicated problem for a 
multiple layered solid half-space. So in this paper we 
present a new counting method to deal with this prob- 
lem. In addition, some new numerical treatments are 
developed: 1) Automatic generation of the travel paths 
and the arrival times of various rays, 2) Automatic gen- 
eration of the products of the reflection coefficients, and 
3) An integration method for head wave rays. 

The conditions for producing various head waves, sur- 
face waves, and interface waves are also examined. 
These conditions are determined from different singu- 
larities in the integrand. 

A FORTRAN program has been developed for nu- 
merical computations of the response for any choice of 
materials for the layer and substrate. The computed re- 
sults for the case of a plexiglass layer and glass substrate 
show that changes of the interface bonding condition 
have a great influence on the behavior of transient waves 
when both the source and detector are located on the top 
surface of the layer. For example, some of the first few 
head waves and regular reflected rays change their po- 
larities and amplitudes when the interface bonding con- 
ditions change. This phenomenon can be used to infer 
the quality of the interface bonding of materials in ultra- 
sonic nondestructive evaluation. In addition, results 
from this fundamental solution are expected to provide 
insight into the study and optimization of probing tools 
such as acoustic microscopes and in applications rang- 
ing from the study of coatings to geo-exploration. 



2. Governing Equations and Boundary 
Conditions 

Consider an elastic structure consisting of an an- 
isotropic homogeneous layer of thickness 2h on a homo- 
geneous half-space as shown in Fig. 1 and suppose there 
is a point force source of step function time dependence 
inside the layer. Then the fields in the layer satisfy the 
equations of motion 



VACUUM 



LAYER 



2h 



(0.0.Z) 



-** X 



LOWER HALF-SPACE 

Fig. 1. Schematic representation of a layered half-space. 
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d<T, 

dx, 



i +f,d(x l )8(x 2 )6(x 3 -z)H(t)- 



h < jc 3 <h 



, d 2 u\ 



(1.1) 



where (0, 0, z) are the position coordinates of the 
source. The origin of the Cartesian coordinates is at the 
center of the top layer; u\, o\j, and p 1 are the displace- 
ment components, the stress components, and the den- 
sity of layer /, respectively; f t , H(t), and 8(x) are the 
components of the point force source, Heaviside step 
function, and Dirac delta function, respectively. The 
summation convention is used. 

The stress and displacement gradient in the layer are 
related by the generalized Hooke's law 



!_ ] d U k 

0-8 - Ciik( Jx~ ( 



(1.2) 



where c l ijU are the elastic constants in the layer. 

Similarly, the field equations and stress and displace- 
ment gradient relation in the lower half- space can be 
written by replacing "I" with "II", thus 



d a,/' „ dW 1 
d Xj d t 



(1.3) 



where a- n ,y, «",, and p" are the stress components, the 
displacement components, and the density in the half- 
space. 

The stress and displacement gradient in the half-space 
are related by 



_ ii _ „n d u k 

<Jij — c ijkt 



dx( ' 



(1.4) 



where c a gkt are the elastic constants in the half-space. 
In addition, the solution should also satisfy the follow- 
ing conditions: 



and 



u k l - u k a - 0, f<0, 
a-i/ = ay 11 = 0, t < 0, 

u k n = 0, x 3 -> -co. 



(1.5) 



(1.6) 



In what follows we consider two cases of the interface 
bonding condition, whereas the top surface boundary 
conditions remain the same. 



Case 1. 

Suppose the interface is "welded", whereas the top 
surface condition is traction free; we have: 



<Tb =0, x 3 = h , 
(To = Ob, x 3 = —h, 

U t l = U; 11 , X 3 = —h. 



(1.7) 
(1.8) 
(1.9) 



Case 2. 



Suppose the interface is intimately "liquid" coupled, 
while the top surface condition is again traction free. We 
have: 



cr f3 = 0, x 3 = h , 

oj 3 = ol\, x 3 = — h , 

oh = CT23 = 0, x 3 = —h, 

o{\ = ol\ = 0, Xi = -h, 

u\ = uf, x 3 = h . 



(1.10) 
(1.11) 
(1.12) 
(1.13) 
(1.14) 



3. Solution Method 

The outline of the solution procedure for the problem 
can be described as follows. First, introduce the Green's 
tensor and take its Fourier transform in time and space; 
then, expand the transformed Green's tensor according 
to the eigenvector of the Christoffel matrix, and decom- 
pose the fields into downgoing waves and upgoing 
waves; third, use boundary conditions to iteratively get 
the solution in the transform domain in a form of 
"generalized ray" series; finally, use the Willis inversion 
technique to get the solution. 

Defining the "Heaviside Green's tensor", G,y, the dis- 
placements in the layer can be expressed as 



ul-Gijfj. 



(1.15) 



Substituting Eqs. (1.2) and (1.15) into Eq. (1.1) gives 



-ijtf 



d G kp 
dXj dX] 



+ 8 lp 8(x 1 )8(x 2 )8(x i -z)H(t) = p 



d Gjp 

dt 2 ' 

(1.16) 



where §,„ is the Kronecker delta. 
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The Green's tensor in the layer may also be expressed 



G = G+ G\ 



(1.17) 



where G™ is the infinite body Heaviside Green's tensor 
and G 1 is the "image" tensor in the layer formed from 
the waves reflected on the boundaries x 3 = ±h. All ma- 
trices are denoted in bold capitals and vectors in bold 
small letters. 

Define matrices K L (co,g) and C L (£) with components 

Kk((o, x) = p L co 2 8 ik - 44^, (1.18) 

ck(x) = chS, (1.19) 

L = I, II, 

where the vector g= (£, £ 2 , &)• 

In terms of Eqs. (1.15)-(1.19), we obtain the equa- 
tions involving G* and G 1 which satisfy 



&[ f t , V )G" = I8(x 1 )8(x 2 )8(x i - z)H(t), (1.20) 



and 



W|,V)G<-0. 



(1.21) 



where / is the identity matrix. The Green's tensor, G", 
in the half-space satisfies 



(1.22) 



^,V)G« = 0. 



And the boundary conditions corresponding to cases 1 
and 2 become: 

Case 1. 

C\V)(G X + G 1 ) = 0, x 3 = h, (1.23) 

C'(V)(G =C + G 1 ) = C u (V)G n , x 3 =-h, (1.24) 

G+G l = G l \ x 3 =-h, (1.25) 

where C'(V) and C n (V) are the matrix operators with 
components 



dXf 



ck<y) = cku^, l = i,ii. 



(1.26) 



Case 2. 

C\V)(G X + G 1 ) = 0, x 3 = /i, (1.27) 

I l C\V)(G<° + G 1 ) = I l C 1 \V)G 1 = 0, x 3 =-h, (1.28) 

^'(VXG^ + G^/zC'^G^O, x 3 =-h, (1.29) 

/sC'tVXG^ + G'WjC"^^ 1 , x 3 =-h, (1.30) 



I 3 (G x + G l )=I 3 G u , x 3 =-h, 



(1.31) 



where 



"1 





0" 











J) 









h = 



"0 





0" 





1 















(1.32) 



(1.33) 



and 



"0 





0" 

















1 



(1.34) 



4. Ray Expansion 

We follow the method developed in Refs. [2,3,4], but 
provide only an outline here. Defining the Fourier trans- 
form of G * by 



G"(co,i)= dn \ dx l dx 2 dx 3 G"(t , x)exp[i(gx — cot)], 



from Eq. (1.20) we then get 



K\w, &6'(w, i) = — /«p(i6z), 

(O 



(2.1) 



(2.2) 



where a> is taken to have negative imaginary parts while 
£ = (£i, £2, 6) nas rea l components. It is easily shown 
that u is an analytic function of co in the lower half of 
the w -plane, and its inverse transform as a function of 
time t is equal to zero when t < 0. 

In order to express the inverse of the matrix K\u>, £) 
in terms of its eigenvectors, we consider the Christoffel 
equation 



K l (a>, i;)u r = A T (u>, i;)u T , 



(2.3) 
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where A(<w, £ ) and u x are the eigenvalues and eigenvec- 
tors. 

For given real <w, £, and £ 2 , the equation 



det £'(<«>, £) = 



(2.4) 



has six roots £ 3 = £ 3 N («, £ a ), iV = ±1, 2, 3, a= 1, 2, 
which may either be real or occur in complex conjugate 
pairs. By means of the concept of Riemann surfaces, the 
six roots may be considered to define a single-valued 
algebraic function &(o>, £«) by Eq. (2.4), if <w is allowed 
to range over the six sheets of its Riemann surface [2], 
It can be shown by analytic continuation that when 
Im(co) < 0, the algebraic function £ 3 has positive imagi- 
nary parts on the three sheets of N= — 1, —2, —3 and 
has negative imaginary parts on the other three sheets of 
/V=+l, +2, +3. 

Now let us consider the eigenvalues. Since A N = 
p\co 2 — (On 2 ), where w N is inverse to £s N (a), £«), let 
£} N (,o)n, £») or 6etK\o) N , £) = 0, and we obtain the six 
roots ±(o N , N = 1,2, 3 and therefore the three eigenval- 
ues A N . 

Normalizing the eigenvectors so that w m w„ T = S mn , we 
have 



I=^u r uJ=UU T 



where the matrix U consists of the three column vectors 
u„ while t/ T is the transpose of U. As a result we have 

3 

G"(u, £) = - -2^"V«r T exp(i6z). (2.6) 

Using symmetry with respect to the x 3 axis, and tak- 
ing the Fourier transform of G 1 from (t,Xi, x 2 , x 3 ) to (w, 
£i, £2, * 3 ) we obtain for the fields in the layer 

G(«, £) = df dxidx 2 G I (f,x)exp[i(£ a x„ - <wf)L 

"" (2.7) 



where 



£ a X„ — £1X1 + £2X2. 



(2.8) 



In what follows we consider the cases of downgoing 
waves and upgoing waves in the layer respectively. 

1. Downgoing waves. 

Consider the case x 3 < z . Here z is the position of the 
point force source. Taking the inverse Fourier transform 
of Eq. (2.6) about £ 3 gives 



2ttw 



G"(w, £, £ 2 , x 3 ) 

3 

d£ 3 2^i~y" r T exp[-i£ 3 (x 3 - z)]. (2.9) 



When x 3 < z , x 3 — z is negative. Then the integral can be 
evaluated by closing the contour in the upper half of the 
£ 3 -plane and by using Cauchy's residue theorem. This 

gives 



G"(w, £, £>, x 3 ) 

_ J_ y « M M M r exp[-i£f~(x 3 - z)] 
-<o£ fdA((o, £) 



5 6 /f 3 =ff 



(2.10) 



where £ 3 M ~ = £ 3 M ~(<w, £, £ 2 ) above, which are located in 
the upper half of the £ 3 -plane when lm(co)<0. The 
subscript or superscript "— " denotes downgoing waves. 

2. Upgoing waves. 

In the case x 3 > z , the contour of integral Eq. (2.9) can 
be closed in the lower half-plane. Similarly we obtain 



(2.5) C + -(», ft, fe, x 3 ) = I i ^exp[-ig(x 3 - z)] 



MQ,£) 
d 6 /f 3 =# 



(2.11) 



where £ 3 A ' + = £ 3 Ar+ ((w, £, §2) are the three roots of Eq. 
(2.4) which are located in the lower half of the £ 3 -plane, 
when Im(io) < 0. The subscript or superscript "+" de- 
notes upgoing waves. 

The general solution of G cc (w, £, £2, X3) should in- 
clude three downgoing plane waves and three upgoing 
plane waves. The reflected waves in the layer should also 
include the downgoing waves and the upgoing waves, 
but there exist only the downgoing waves in the lower 
half-space. Then we respectively obtain 

G°°(<u, £, &, x 3 ) = G-(u>, £, £ 2 , Xi) + G+(w, £, £ 2 , X3), 

(2.12) 

3 
G' = 2 M MV M T exp[-i£ 3 M_ (x 3 - h)] 

3 

+ ^u N v N expt-i^fe - h)], -h<x, <h, (2.13) 

3 

G" = 2^ w P T exp[-i£T(x 3 - h)], x 3 < -h, (2.14) 
where qp represents the eigenvectors of an infinite 
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body in the lower half-space while £ p ~ 3 represent the 
roots of Eq. (2.4) when the superscript "I" is replaced by 
"II", for the same reasoning as mentioned above, which 
roots locate on the upper half-plane of £s-plane when 
Im(co) < 0, and v~ M , v + N , and w~ P are the coefficients to 
be determined. 

Further matrix notation is introduced as follows, 



U = \u\ , u 2 , u 3 ], 

v~ = [vr, V2, vj], 



(2.15) 
(2.16) 



with corresponding definitions when superscript " — " is 
replaced by "+" and 



Q -[qi,qi, q-i ], 


(2.17) 


N~ = [wi, w 2 , wj], 


(2.18) 



Then the previous equations can be written in more 
compact form as 

G a ((o, & &, xi) = H(z - x,)U-p-(z -Xi- h)D-U- J 

(2.19) 
+ H(x 3 ~ z)U + P + (z - x, + h)D + U +T , 

where P , P + , D , and D + are diagonal matrices with 
components 



(p-)((xd = 5/exp[i^-fe + h)], 



where 8[ is the Kronicker delta. Also, 



(2.20) 



(P + )((*i) 


= S/exp[i£f + (x 3 


-h)], 


(2.21) 


<p-)l=8t 


[^l 


J' 


(2.22) 


{D + )i=8[ 


H^L 


-&) ' 


(2.23) 



and 

G 1 = U-p-(-x 3 )V- J + U + P + (-x 3 )V +T , (2.24) 

G" = Q-S-(-x 3 )W- T , (2.25) 

where S~ is the diagonal matrix with components 

GT)((x 3 ) = S/exp[i^-(x 3 - *)]. (2.26) 



Case 1. 

If we are only interested in the fields in the layer, 
substituting Eqs. (2.19), (2.24), and (2.25) into the 
Fourier transformed equations corresponding to the 
boundary conditions Eqs. (1.23)-(1.25) and eliminating 
the coefficient matrix W, we get 



B\P + V +1 = - B\P\z)D + U + 



(2.27) 



and 



B + 2 V +1 - B^PV 1 = - B2P(z)DU T (2.28) 

where 

B\ = C\g)ir, (2.29) 

B + 2 = [C\i) - C n {m'{Q~)~ l W\ (2.30) 

BT = C\£)ir, (2.31) 

B- 2 = [C"(f> - C n (£)Q-(Q-y ] ]U-, (2.32) 

P + =-P + (-h), (2.33) 

P = -p-(h). (2.34) 

As h approaches infinity the exponential functions P + 
and P on the left-hand sides of Eqs. (2.27) and (2.28) 
approach zero because the exponential factors have neg- 
ative and positive imaginary parts, respectively. There- 
fore, Eqs. (2.27) and (2.28) uncouple into two indepen- 
dent equations and the problem becomes two separate 
problems for half-spaces x 3 <h and x 3 > — h. In addi- 
tion, the larger Im(co) becomes the smaller exponential 
functions P* and P , and so the equations can be solved 
iteratively to give a uniformly convergent series of 
"generalized rays" for given £i, t- 2 , and the real part of a>. 
Finding V~ T and V +T from Eqs. (2.27) and (2.28) by 
means of an iteration technique similar to Ref. [5] and 
substituting the results obtained into the previous rele- 
vant equations we get the general field in the layer 



G = G + G + , 



(2.35) 



where 



G = H(z - x,)U-p-(z -Xi- h)S- 
U'P'(-x i )R + P + (z)S + - irp-(- Xi )R + P + R-p-(z)S- 
- U-p-(-x 3 )R + P + R-p-R + P + (z)S + - ... (2.36) 
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and 



G + = H( Xi - z)U + P + (z - xi + h)S+ 
U + P + (-Xi)R-p-(z)S- - U + P + (-Xi)Rp-R + P + (z)S + 



U + P + ( -Xi)R-pR + P + Rp- (z )S~ 



(2.37) 



which are the downgoing waves and the upgoing waves 
in the layer, respectively. Their physical meanings are 
obvious. If we notice that R + and R are respectively the 
matrices of the reflection coefficients at the top surface 
and the lower surface of the layer, and the diagonal 
matrices P + and P represent phase delays, then the first 
term of Eq. (2.36) represents the downgoing direct rays 
from the source to the observation point (detector), the 
second term describes the rays which reflect once at the 
top surface of the layer and so on. Similarly, we can 
explain each term of Eq. (2.37). Using the concept of an 
infinite linear array of image sources, then, with the 
exception of the first (direct wave) term, all the other 
terms of Eqs. (2.36) or (2.37) can be considered to be 
produced from the corresponding image sources located 
above or below the layer. The definitions of the notation 
of Eqs. (2.36) and (2.37) are as follows: 



If = (BIT 1 , L + 2 = (B + 2 y\ 

R = L2D2 , R = L\ 2>i, 

S~=D~ir T , S + = D + U +T . 



(2.38) 
(2.39) 
(2.40) 



A typical term of Eq. (2.36) except for the first term 
may be written as 



U-p-(-x 3 )R + P + R-p-...R ( - ir P ( - l, (z)S ( - l) , (2.41) 

where k counts the number of P 's, R 's or 5"s in the term. 
Interchanging superscripts "+" and "— " in Eq. (2.41) 
gives the typical term of Eq. (2.37). 

Suppose an element of the matrix of a typical term of 
Eq. (2.41) is written in the form 

F(a>, &xp[-K€aXa -<»t + A,£ M + A 2 $)]. 

Then the inverse Fourier transform of Eq. (2.36) gives 

(G-) m ' ^3 \ d*>J"Jd6d62>(«. 

exp[-i(f«jc« -o>t+ A^ 3 M + A 2 g)], (2.42) 

where (G~) m( is the element of the matrix G and F(w, 
g) is a homogeneous function of degree —2, because the 



elements of the matrices U ± and R ± are homogeneous 
functions of degree zero, while the elements of the ma- 
trices S ± are homogeneous functions of degree —2. In 
addition, since the series of "generalized rays" is uni- 
formly convergent, we can invert the integration term by 
term using the Willis inversion method (Appendix A) to 
obtain 



(G-r = =± 1 I ds rr^fW . (2-43) 
4tt J -t + Ai£ 3]fi + A 2 £ 3 ,fl 

1/11=1 

...to 



where (G ) mt is the element of G , while the subscript 
fl denotes d/dfl and 



O-— -is. 



(2.44) 



fl in the integrand is now taken as the root in the lower 
half-plane of the equation 

-flt+ r, a x a + X^ifl, i, a ) + kigiQ, Va) = ei, (2.45) 

where e is an arbitrary infinitesimal number. 
For the elements of G + , we have similar results. 

Case 2. 

In this case, substituting Eqs. (2.19), (2.24), and 
(2.25) into the Fourier transformed equations corre- 
sponding to the boundary conditions Eqs. (1.27)-(1.29), 
we get the same Eqs. (2.27) and (2.28) when B + 2 and B 2 
are replaced by B* + 2 and B*~ 2 , respectively, and 

/?•* = {cx& - {hc\m~) 

[{h+I 2 )C l \^)Q-+hQ-V}U ± . (2.46) 



Defining 



L 2 + = (fl 2 *T\ 

R = L<2 $2 1 



(2.47) 
(2.48) 



and using the matrices with "*" instead of the corre- 
sponding matrices without "*" in Eqs. (2.35)-(2.36), we 
immediately obtain the solution of case 2. The inversion 
of the transform is the same. 

Up to now we have not considered the properties of 
the materials, so all the previous results can be applied 
to the general case of arbitrary anisotropic materials. 
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5. Special Case of Isotropic Materials 

When the layer and the lower half-space consist of 
two different isotropic materials, then the components 
of K L (co, £) become 

= (frco 2 - 1X^)8* - (A L + ^ L )t&, L = I, II. (3.1) 

where A L and p L are Lame's elastic constants of the 
materials in the layer (L = I) and the lower half- space 
(L = II). 

The roots of the equations 



:irc 



det# L (w,£) = 0, L = I, II 

£f" L = ±p L (once), 
g cL =±q L (twice), 



where 






or 
at 



cl 



L = I, II, 



L = I, II, 



at = (A L + /j, l )/pl, L = I, II, 
ct = I^JPl, L = I, II, 



e=e+&. 



(3.2) 

(3.3) 
(3.4) 

(3.5) 

(3.6) 

(3.7) 
(3.8) 
(3.9) 



Making a comparison with the previous sections gives 

£ 3 M - = p h q h qi 

6 = -pi, -qi, -qi. 



(3.10) 



Three roots, p h q h and q h of the six roots in Eq. (3.10) 
are in the upper half of the £-plane, whereas the other 
three roots, —p h —q h and — q u are in the lower half of 
the f 3 -plane. For the details of the distribution of these 
roots, refer to Ref. [2] 

The normalized eigenvectors corresponding to the 
roots of Eq. (3.10) are 



6 

+pi. 



(e + (.Pifr 



(3.11) 



u 2 = 



~-&qi 

-e. 



(e+(qi) 2 r m r, 



(3.12) 



and 



"3 



-fer 1 





(3.13) 



For the lower half-space there exist only the downgoing 
waves; the following three roots should be chosen as 



(3.14) 





& p - = 


Pa, qa, qa. 


The corresponding eigenvectors are 


<?r = 


"6" 
6 

-Pn- 


ie + (Pafr m , 


qi = 


&qn 

-el 


[e + (qa) 2 r" 2 ^, 



(3.15) 



(3.16) 



and 



«3 = 



-fer 1 





(3.17) 



The matrix operators C'(£) and C n (£) should always 
operate on the respective regions of eigenvectors first to 
obtain the proper stress components. We have 



: 2/xi/?igi frig? - g% .jMqS' 



cx&u* 



Vpjve VqTVe-t z 

+ 2/j,iPig 2 iijjqt ~ g 2 )fc .jMq&_ 



Viiqt ~ j 2 ) , 2/Ai£if 



f 



_ Vpj^e Vqjve-t 



(3.18) 



c\€)u ± = 



: 2ju, n pngi [Anjqti ~ g 2 )gi , jAaqa& 



VpJVf VqJVf-Z € 

+ 2/Aiipiigz p-niqh ~ g 2 )fc ^ Mngngi 

Vpjve Vqt +e-z ~ f 

Mn(gn ~ g 2 ) . 2p, n qn^ 2 



VpITIF V^ + ^-f ' _ 

(3.19) 
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Substituting the relevant matrices into Eqs. (2.39) and 
(2.48) we obtain the reflection coefficient matrices R + , 
R , and R*~ respectively. 

R + = UVl (3-20) 

R + n = R + 22 = [(q? ~ ff - Mq.eVKq, 2 - £ 2 ) 2 + 4 Ml £ 2 ], 

r + u = [4q&q? - e)(P? + ey i2 v{(qi 2 + em* 2 - e? 

r + 21 = - [4 Pl t(q? - exqi + er-viipi + mc* 2 - e? 
+4 W ,a), 

^?33 — — 1, 

Rn — Rii = R31 = ^32 — 0, 

R = [R$], (3.21) 

Ru = [-A$Ai + AJAQ/An, 

Rn = [-At At + A\AT}IA m 

R 21 = [AtA 2 + AZAjyAn, 

R 22 = [At A; - 4J4r]/4 n , 

R33 = (Mi<?i - ^iiqn)l{^iqi + (inqn), 

R13 — R21 = R31 = R32 = 0, 

where 

Af = { ±2 f ji l qi + Uniqiiql + f) 

- qi (2 P2 q 2 + e- qZWtowi + £ z ))/(q 2 + £ 2 )" 2 , 

Af = { +2yu,i/?i + ixii[-p 2 (qi + i 2 ) 

± Pl (2p 2 q 2 + e~ q 2 W(p 2 q 2 + f*)}/(p? + ?) m > 

Af = {fliiqi - £ 2 ) + p, n [(2p 2 q 2 + £ z - q}) 

± piqfal + IWbw + e)}/(j> 2 + ey 2 , 

Af = {friqi - f) + IX„[(±p 2 qi ~ l)(qi + 1) 

+ (2p 2 q 2 + 2£ - qi)]/(p 2 q 2 + f 2 )}/^! + £ 2 )" 2 , 

A n = At At ~ At At, 



and 



R'~ = [Rv~], (3.22) 

Rn = [-4 1 * + 4r + 4r4 4 * + ]/4n, 
Rn = (-4; + 4r + 4 4 * + 4D/4n, 
Rn = (A; + At + 4 2 * + 4 3 *-)/4,* I , 

Rl 2 = [a; + a;- - 4 2 * + 4 1 *-]/4 I ;, 

R33 = 0*i?i ~~ ^nqn)l{^iqi + jttii^n), 
Rn = R 2 3 — R31 = ^32 =0, 



where 



4/^±2^ 1 + ^-^^ W + l,^ 



A? = 'Ip^lipl + l) 1 

4., - 1 M* 2 - 1) ± ™ i[iq :rJ f + 4p2q2] W; + 1 <" 



Ar = p, 2 



p 2 {q 2 + \) 
4g 2 (p 2 - g 2 ) p 2 u> 2 



p 2 (q 2 + 1) p 2 



In what follows let us first consider the interface con- 
dition for case 1, i.e., the welded interface. In the case of 
isotropic materials the inverse transform Eq. (2.43) 
takes the following form: 



(fry* = 2(- 1)*[1 + 8l{H{z - X3) - 1)] A 

Sty {ii,V> 

en 



8=1.2.3 J 
j=1.2,3 1,1=1 
d=ll 



where j and g take 1, 2, and 3 and respectively represent 
the type of the final trip and the initial trip of the ray 
path to be a P ray, SV ray or SH ray in the layer, and 

4>(0, 77) = —flt+ rj a x a — 2hkipi — 2hk 2 qi, (3.30) 

where k~\ and k~ 2 are the numbers of times the layer is 
traversed by a P ray and S ray, respectively. They may 
take fractional values when the source or the detector is 
not on the top surface or the interface. Also, 



{D-)i=-(2p l nb k )- l S[, 



(3.31) 
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where 



b\ = alpi, b 2 = b 3 = c?q h 



(3.32) 



R (j , g , k) represents the sum of the products of all the 
reflection coefficients produced by the rays which travel 
from the source to the observation point (detector) and 
touch the top surface and the bottom surface of the layer. 
These rays have the same final trip j and the same initial 
trip g as well as the same total number of layer traverses 
k of the layer as shown in Fig. 2. I.e., they have the same 
arrival time. The concrete expression of R~(j, g, k) will 
be discussed in the next section. 



DETECTOR 



SOURCE 




Fig. 2. A typical ray path within the layer; p and s respectively 
represent the sections traveled with p wave speed and s wave speed. 
i? ps , i? pp , etc., respectively, represent the reflection coefficients of the 
top surface and the lower surface of the layer; the superscripts or 
subscripts represent the wave mode before and after reflection. 

Similarly, we have 



(GT'=i (-!)*{! + 8[[H(x 3 - z) ~ 1]} -k 



v i (u+)rR + (j,g,k)[(u^n T mD ( - 1 




7=1,2,3 |,| = 1 — 


(3.33) 



where 

4> + (fl, tj) = -fit + T) a x a - IhkXpi ~ 2hktq h (3.34) 

{D + Y k =-(2 Pl nb k y x 8 t k . (3.35) 

The meanings of the various quantities in Eqs. (3.33)- 
(3.35) are the same as before or are similar to that of the 
quantities in Eqs. (3.29)-(3.32). 
Using 

<£(A V) = 0, (3.36) 

then Eqs. (3.30) and (3.34) can be written as 

2Mf 2hk 2 ^ 



c) 4> ± \ _! 

d ii | ^= fi 



'fa-^a 



Pi 



<7i 



5.1 Detector on the Top Surface 

Usually, the detector (observation point) is put on the 
top surface of the layer; then Eqs. (3.29) and (3.33) can 
be further combined. If the detector was considered to 
be located an infinitesimal amount below the upper sur- 
face, we should have to take into account at almost the 
same instant the wave going upward before reflection 
and the reflected wave going downward which come 
from the same ray and have almost the same phase 
function. Substituting into Eqs. (3.29) and (3.33) with 
k = n + 1 and k = n + 2, respectively, and defining 



r=u + - ur\ 



(3.38) 



and substituting the previous relevant expressions into 
Eq. (3.38) we obtain 



Here 



(3.37) 



r= [y,jl (3.39) 

yu = lvApiqi(qi+l)VlPi-+l] l,2 A h 

jn = [~Vi2qi(qi ~ 1)(<?i 2 + 1) 1,2 ]M„ 

yn = -2172, 

72i = Im^pMqi + i)M/>. 2 + i] 1/2 4, 

722 = [~V22qi(qi ~ l)(qi + 1) 1,2 ]M„ 

723 = 2l7i, 

T3i = [-2 Pl (qi ~ l)(?i 2 + WlPi + 1]" 2 4, 
732 = [-4p I ^, 2 +l) 1/2 ]/Z\„ 

733 = 0, 

A l = {qf - l) 2 + 4 Pl q h 

Then Eqs. (3.29) and (3.33) can be combined and writ- 
ten as one formula: 

„_n 8=1,2.3 . 

j=l,2.3 |,| = | 
0=0 

(-IfyfRU, g, n){U^y g D4>{k u k 2 )P(g)ds, (3.40) 
where 

4> = —ry + d — kipi — k 2 qi, (3.41) 
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T=Cit/2h, y = Qlci £; a = x a /2h, d=^ a rj a , 
Pl = (a 2 y 2 ~ l)" 2 , qi = (y 2 ~ l)" 2 , g=z/2h, 
ai-cja h k p = n p + 8f\- - g), j8=l,2. 



D<p(k u k 2 ) ■ 



1 



16tt 2 p\h\ d + — H — 
Pi qi 






(3.42) 
(3.43) 

(3.44) 



5.2 Source and Detector on the Top Surface 

Similarly, when the source is also located on the top 
surface of the layer then Eq. (3.40) may be written as 

„_n *=1,2,3 J 

n ~ V j=l,2,3 \ v \=\ 

(-lTyrR(j,g,n)(iliiD ( t>(k 1 ,k2)P(g)ds, (3.45) 
where 

4>= —Ty + d— kipi — k 2 q t . (3.46) 

(if/)' g is a component of the following matrix, 

y = R + (U + f- (£T) T . (3.47) 

Obviously, substituting the corresponding quantities 
with "*" into Eqs. (3.29), (3.33), (3.40), and (3.45) gives 
the formulas corresponding to case 2, i.e., the liquid 
coupled interface. 

If we are only interested in the early time or short 
time arrival of the signal received by the detector, only 
the finite terms of all the previous integrals need be 
computed, because the reflected rays with more reflec- 
tions do not have enough time to arrive at the detector. 

Following Ref. [3], we introduce the following trans- 
form 



171 = cos0 = dlx, 

V2 = & md= ±{x 2 - d 2 ) m lx 

ds = dd 

y{6)=y(-B) 



(3.48) 



where x = \x\ and 0is the angle measured around 1 77 1 = 1 
from the direction of the vector x . Note that if y is a root 
of 

-Ty + d- ka(l - a 2 y 2 ) [ ' 2 - k 2 i(l - y 2 )" 2 = 0, (3.49) 



for a given value of d, then — y will be a root of the 
following equation when d is replaced by —d: 

ry-d- kii{\ - a 2 y 2 ) m - fei(l - y 2 ) m = 0. (3.50) 

Thus, if d and y form a solution to Eq. (3.49) or (3.50), 
then d, and y (they are the complex conjugates of d, y) 
are solutions to Eqs. (3.50) or (3.49), while — d, and — y 
are solutions to Eqs. (3.49) or (3.50). When all the 
branch cuts for the square root functions in the above 
equations were taken along the real axis between each 
pair of branch points, using the above relations Eqs. 
(3.48)-(3.50) and, finally, the alternate form cf>= si for 
the integration locus, then every term of the integrals 
Eqs. (3.40) and (3.45) can be simplified to the form 



I[y(d),d]-^ 

-x-i0 \ X 



dd 



d 2 ) 1 



(3.51) 



where the integrand I[y(d), d] represents the integrand 
of each term, including all the factors of Eq. (3.40) or 
(3.45). The denominator (x 2 — d 2 ) 1 ' 2 and the constant 
factor 2 are introduced by the variable transform Eq. 
(3.48) and the change of the integration path. All the 
singularities of the integrand of Eq. (3.51) are located in 
the upper half or on the real axis of the d -plane. For a 
detailed derivation of Eq. (3.51), refer to Ref. [3]. 



6. Automatic Generation of Ray Paths 
and of Corresponding Products of 
Reflection Coefficients 

In general, there are many rays which arrive at the 
detector at the same time owing to the multiple reflec- 
tion and the mode conversion of the incident P ray or S V 
ray emitted from the force source in the layer, and there- 
fore it is not necessary to separately compute each term 
in the integrals mentioned above. The question of how 
many rays will arrive at the detector at the same time for 
a given configuration is a problem of combinatorics. 
Although explicit formulae for summing up the rays 
with the same arrival time to form a single integral can 
be derived for the case of a plate because the reflection 
coefficients at the top surface and bottom surface are the 
same [3,12], similar formulae are not easily derived for 
the present case of a layer on a substrate. Even though 
these formulae can be obtained, the actual program 
would also be rather complicated to carry out for a 
summation with reflection coefficients at the top surface 
which are different from those at the lower surface. In 
the following we will give a special counting scheme for 
this problem. First, let us change the notation in Fig. 2 
and we obtain Fig. 3. Notation "1" and "0" in Fig. 3 
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DETECTOR 



R(0,1,2) 



R(1,0,2) 




SOURCE 



R(1,0,1) R(0,0,1) R(1,1,1) 



Fig. 3. A typical ray path within the layer; p and s respectively are 
replaced by and 1; fi ps , etc., by i?(0,l,2), etc. 

respectively corresponds to notation "s" and "p" in Fig. 
2. As a result it is not difficult to understand the mean- 
ing of R(0,1,1), R(l,l,2), etc. The first argument de- 
notes the type of incident ray while the second one the 
reflected ray. The third arguments, "1" and "2", respec- 
tively, represent the lower surface and the top surface of 
the layer. Obviously, the configuration of the rays in Fig. 
3 may be expressed in a sequence consisting of "0" and 
"1" as shown in Fig. 4. On the other hand, the sequence 
in Fig. 4 may also be considered as a binary notation of 
some integer. 







1 











Fig. 4. An example of a binary sequence, that uniquely characterizes 
the ray path in Fig. 3. 

To find how many rays travel through the thickness of 
the layer three times with p-wave speed and three times 
with s-wave speed, therefore, is the same problem as to 
find how many different six-bit binary integers can be 
constructed by three l's and three 0's. Using a bit-test 
program it is easy to get the count and find all those 
binary numbers that correspond to all the ray paths with 
the proper p-wave and s-wave sequences. An outline of 
the procedures for assembling the proper sequence of 
the reflection coefficients for each ray path is given in 
the following for the case of the source and the detector 
located on the top surface of the layer. 

1 . For given k x and k 2 , the number of layer traverses 
with p-wave speed and s-wave speed respectively, there 
is in general more than one ray path which will have the 
same arrival time. But each ray may have a different 
reflection sequence. From Figs. 3 and 4, the sequence 
for each ray path has a corresponding N-bit binary num- 
ber representation of k\ l's and k 2 0's; N= ki + k 2 . 

2. Using the bit-test program, we can determine IP, 
the number of distinct binary sequence with k t 1 's and 
k 2 0's, store all such N-bit binary numbers in an array, 
IA, of dimension IP. Each binary number will have k t l's 
and k 2 0's, and they are all different. 

3. For each binary number in IA we can construct a 
product of the reflection coefficients according to Fig. 
3. 



4. Summing all the products thus obtained gives all 
the rays which arrive at the detector at the same time. 

The procedure outlined above is rather efficient in 
terms of computation speed as long as N is limited to a 
small number, say less than 16. The results in the case 
of a plate without the lower half space have been 
checked against the results computed using explicit for- 
mulae previously derived [12]. These same formulae for 
a plate cannot be used for the case of a layer on a 
half- space. 



7. Singularities and Wave Front Arrivals 

The integrand of Eqs. (3.40) or (3.45) may contain 
singularities, and they can be divided into three cate- 
gories that respectively correspond to three different 
types of wave front arrivals: 1) The body waves relevant 
to the regular reflected rays which are determined from 
the singularities of dcf>/dy = 0, 2) The interface waves, 
including Rayleigh surface waves, Stoneley interface 
waves, and the other possible leaky waves, which are 
determined from the singularities of A\ = or A B = in 
the denominators of the reflection coefficients, and 3) 
The head waves determined from the singularities of the 
branch points of the square root functions. 

Following the analysis in Ref. [3] we will show that 
all the above mentioned singularities of the integrands in 
the y -plane will be located in the upper half or on the 
real axis of the y -plane and their mapping onto the 
<i-plane, which is determined from Eq. (3.48), will lo- 
cate them in the upper half or on the real axis of the 
<i-plane. So the integration path will never touch the 
singularities of the d -plane when we use formula (3.51) 
to compute the Green's functions and choose the inte- 
gration path slightly below the real axis. This is indeed 
a main advantage of the Willis method. 

Let us first consider the head wave arrivals which are 
related to the singularities of the branch points. 



8. Head Waves 

Case 1. Welded interface. 

In this case we have the following four square root 
functions: 



p^iaty'-rr, qi = (a 2 y-iy\ 



(4.1) 



Pll = (a 3 y " 1) , <Zn = {ccly 1 - l) 1 ", 



where 
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«i = Ci/ai, a 2 =\, a 3 = Ci/an, a 4 = c l /c u . (4.2) 

Their branch points are respectively 

yi =±ai\ y 2 =±l, yi = ±a,\ y 4 = ±a 4 \ (4.3) 

For convenience of discussion, without losing gener- 
ality we assume the source and detector are located on 
the top surface of the layer and 

an > c u > ci\ > Ci. (4.4) 

Then from Eq. (4.2) we have 

a 3 < a 4 < ai < 1. (4.5) 

The relation between y and d satisfies the equation 

d = ry + /t,i(l - a 2 y 2 ) m + k 2 i{\ - y 2 ) m , (4.6) 

where 

T=tcJ2h, d = x a r] a /2h. (4.7) 

A simple analysis of Eq. (4.6) shows that when all the 
branch cuts were taken on the real axis of the y -plane 
between each pair of branch points, the mapping of the 
branch cuts in the <i -plane will be located in the upper 
half and on the real axis of the d -plane. 

The following gives some examples of head wave 
arrivals corresponding to the branch points. 

1. The head wave arrivals corresponding to branch 
point y t . 

From Eq. (4.3) we have 



and 



yi = \la x = cii/ci, 



y = csct 



(4.8) 



(4.9) 



where 6 is the incident angle. Substituting Eq. (4.9) into 
Eq. (4.8) and writing 8 as d a i we obtain 



sin0„i = Ci/oi. 



(4.10) 



Equation (4.10) is the critical condition satisfied by the 
critical angle B a i which leads to the generation of one 
type of head wave. The head waves of this type is named 
SP*S head waves. Here S and P respectively denote an 
SV ray and P ray in the layer. For a detailed physical 
mechanisms of the generation of head waves, refer to 
Refs. [6] and [7]. 



Substituting Eq. (4.8) into Eq. (4.6) and letting k x = 
we obtain the normalized arrival time of the SP*S head 



J Cl , 7 Cl 
t= d — v k 2 



U\ 



Cl\ 



(4.11) 



For example, when k 2 = 2 it is just the arrival time of the 
SP*S head wave, where SP*S also represents the propa- 
gation path of the head wave ray. It consists of three 
sections, the first section of the notation, S, represents a 
critically incident S ray in the layer, the second, P*, 
represents a P ray propagating along the interface but on 
the side of the layer, and the final S is a critically 
reflected S ray in the layer. Later on, the lower case letter 
p or s with a star will represent a P ray or S ray propa- 
gating along the interface on the side of the half-space. 

2. The wave arrival corresponding to branch point y 2 . 

From Eq. (4.3) we have 

y 2 =\. (4.12) 

The corresponding critical angle satisfies 

sin0 a2 = 1, or 6 a2 = tt/2; (4.13) 

this condition seems as if the S ray emitted from the 
source propagates along the top surface of the layer. We 
denote this ray by S*. It is possible only if k\ = and 
k 2 = 0. Letting y 2 = 1 in Eq. (4.6), we obtain the normal- 
ized arrival time of this ray 



T=d. 



(4.14) 



3. The head wave arrivals corresponding to branch 
point y 3 . 

From Eq. (4.3) we have 



y 3 = On/ci. 
The corresponding critical angle is 

6^ = sin _1 (ci/fl n ). 



(4.15) 



(4.16) 



In this case we may have two kinds of head waves that 
are excited by the same critically incident S ray but have 
different critically reflected S and P rays. With k x = 
and using Eq. (4.15) in Eq. (4.6) we obtain the normal- 
ized arrival time of one kind of head wave 



d v k 2 — 

flu an 



(4.17) 
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When k 2 = 2 we obtain the arrival time of the Sp*S head 
wave. 

Using Eq. (4.15) in Eq. (4.6) we obtain the normal- 
ized arrival time of the other kind of head wave 



t = a h k\ 



flu 



flu 



+ fe 



_£i_ 
' An 



(4.18) 



When ki = 1 and k 2 = 1, it is the arrival time of the Sp*P 
head wave. 

4. The head wave arrivals corresponding to branch 
point y 4 . 

From Eq. (4.3) we have 



y 4 = di/ci. 
The corresponding critical angle is 
a4 = shr'tci/cii). 



(4.19) 



(4.20) 



In this case we also have two kinds of head waves. With 
ki = and using Eq. (4.19) in Eq. (4.6) we obtain the 
normalized arrival time of one kind of head wave 



T = fl h k 2 



flu 



Cn 



(4.21) 



When k 2 = 2 it is the arrival time of the Ss*S head wave. 
Using Eq. (4.19) in Eq. (4.6), we obtain the normal- 
ized arrival time of the other kind of head wave 



t = a — v k< — 



Cn 



Cn 



+ kn 



Ci 



(4.22) 



When ki = I and k 2 = 1, it is the arrival time of the Ss*P 
head wave. 

It may be seen that the branch points y u y 3 , and y 4 
correspond to the arrivals of various head wave rays, 
while the branch point y 2 seems to correspond to the 
arrival of the S ray propagating along the top surface of 
the layer. These head waves all are excited by the S ray 
emitted from the source under the conditions of Eq. 
(4.4). A picture of the wave fronts of the head waves in 
the layer generated by reflection and refraction of a 
spherical P wave and a spherical S wave impinging on 
the interface is shown in Fig. 5. 

In order to describe the arrivals of various head wave 
rays and the P ray propagating along the top surface of 
the layer excited by the P ray emitted from the source, 
we need to change the form of the relations mentioned 
above. We write 



POINT FORCE 




Fig. 5. Schematic representation of the head waves in the layer gener- 
ated by reflection and refraction of a spherical P wave and a spherical 
S wave impinging on the interface. 



_ Q _ a.\fl _ ciifl _ ci\ 

C\ Cl\C\ C\Cl\ C\ 



where 



y p = illai. 
Substituting Eq. (4.23) into Eq. (4.1) gives 



Pl = (tfyl ~ D 1/2 , fli = fe 2 " D 1 ' 2 , 



(4.23) 



(4.24) 



(4.25) 



Pll = (tfyl - 1) 1/2 , a„ = (/3 4 2 y 2 - l) 1 ' 2 , 
where 

j8i = 1, /3 2 = fli/ci, /3 3 = ai/an, /3 4 = ai/c n . (4.26) 
Substituting Eq. (4.23) into Eq. (4.6) we obtain 

fli 



d=r^y v + k l i(\ -v p 2 )" 2 + fci 

C\ 



1 



C\ 



)-'p 



■ (4.27) 



Obviously, the transform of Eq. (4.23) maps the orig- 
inal y -plane onto the v p -plane, whereas the original 
branch points of the y -plane become those of the y p - 
plane as follows: 

y p , = ±l, y P 2=±fc\ yp3 = ±/83 _1 , y v *=±fc\ 

(4.28) 

Similarly, we may discuss the arrivals corresponding to 
the branch points of Eq. (4.28). 

5. The head wave arrivals corresponding to the SH 
ray. 

Since the SH rays do not produce mode conversion 
and their speed is equal to that of the SV rays, the 
arrivals of the head waves excited by the incident SH 
rays emitted from the source are the same as those of the 
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Ss*S head waves. These head waves may be called the 
Hh*H head waves. 

Case 2. Liquid coupled interface. 

In this case the number and the distribution of the 
branch points are the same as those in the case of the 
welded interface. The arrival time of various head waves 
and surface P and S rays are the same as those of the 
welded interface. (The polarities and amplitudes have 
important diagnostic differences; see Computed Results 
and Discussion.) 



9. Interface Waves 

The interface waves include the Stoneley waves and 
the leaky waves which propagate along the interface of 
two different solid media [8]. When one of the media is 
vacuum we obtain the Rayleigh waves propagating along 
the free surface of the solid. The Stoneley waves do not 
attenuate and the leaky waves attenuate when propagat- 
ing along the interface. Unlike the Rayleigh waves, the 
Stoneley waves can exist only if some parameters of the 
media satisfy certain conditions [9]. The singularities 
corresponding to the Stoneley waves come from the real 
roots of the so called Stoneley equation, i.e., A a = on 
the top Riemman sheet of the y -plane, while the singu- 
larities corresponding to the leaky waves come from the 
complex roots of the Stoneley equation on the other 
Riemman sheets of the y -plane and therefore have some 
attenuation [10]. 

Suppose the source and the detector are located on the 
top surface; in this case it can be shown from Eq. (3.49) 
that the mapping of the real root singularities corre- 
sponding to the Stoneley waves onto the d -plane will 
locate them in the upper half of the d -plane, while the 
mapping of the complex root singularities correspond- 
ing to the leaky waves onto the d -plane will never occur 
on the top Riemman sheet of the d -plane. Thus only the 
roots of the Rayleigh equation Aj = will be considered. 
Letting the denominator of the reflection coefficients of 
the top surface be equal to zero, we obtain the following 
Rayleigh equation: 



(qf - l) 2 + 4p iqi = 0. 



(4.29) 



Using the principle of the argument, it can be proved 
that Eq. (4.29) has only two real roots, +y R and — y R . 
Obviously, only the positive real root is of interest. The 
speed of the Rayleigh waves is less than that of the S 
waves, and so y R < 1 . 

Substituting y R into Eq. (4.6) and letting k x = and 
k 2 = 0, because the source and the detector are located 



on the top surface, we have 
d = ry R . 



(4.30) 



Since t and y R are real numbers, the d corresponding to 
the singularity y R should be located on the real axis of 
the rf -plane. The arrival time of the Rayleigh waves is 
determined from Eq. (4.30) if d and y R are known. 



10. Regular Reflected Rays 

Analysis shows that the singularities of d<[>/dy = 
which relate to the regular reflected rays are located on 
the real axis of the y -plane and so the mapping of these 
singularities onto the d -plane will locate them on the 
real axis of the d-plane from Eq. (3.49). It is not diffi- 
cult to understand this consequence if we note that 
y = cscf?, where is the incident angle of the ray, and d 
is the wave front distance satisfying Snell's law. When 
d, a u ku and k 2 are given, we have a group of specified 
rays that have the same arrival time t satisfying Eq. (4.6) 
and the following equation: 



d = 



(a.V 



D 1 " (y l - it 



(4.31) 



Equation (4.31) may be obtained using Snell's law. Ob- 
viously, we must first know the value of y in order to get 
the arrival time r. It is seen from Eqs. (4.6) and (4.31) 
that this is a problem of solving a nonlinear algebraic 
equation in y . Following Ref. [12], the solution for y is 
obtained using a simple, intuitive, iterative technique. 
The arrival time t is determined by substituting the 
solution obtained into Eq. (4.6). 



11. Integration Technique and 
Computation Procedures 

The previous analysis shows that when we use Eq. 
(3.51) and choose the path of integration below the real 
axis of the d -plane, it will not touch the singularities of 
the integrand, which include Rayleigh poles and all the 
branch points of the square root functions. When all the 
branch cuts are taken along the real axis of the y -plane, 
then their mapping in the d -plane will locate them in the 
upper half-plane and on the real axis. This is indeed a 
main advantage of the Willis method. 

The actual procedure for numerically computing the 
Green's tensor can be summarized as follows: 

1. For given material properties and test configuration 
represented by x, the distance between source and detec- 
tor, a time of arrival table is computed first. The arrivals 
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include all the regular reflected rays and all possible 
head waves and Rayleigh wave; each has an associated 
pair of number k x and k 2 , denoting the group of the ray 
paths. The arrivals are sorted according to successive 
time sequences. 

2. To compute a particular component of the Green's 
tensor for a specific time, t, the number of possible 
arrivals, N, can be determined by a comparison with the 
time of arrival table computed in step 1 . N is also the 
number of terms in Eq. (3.45) that need to be computed; 
each term consists of one definite integral. 

3. For each integral to be numerically computed, the 
integrand consists of IP terms and each is a product of 
a unique sequence of reflection coefficients and some 
other factors. Both IP and the sequence series are com- 
puted by calling a bit test program as explained in Sec. 
6. 

4. The numerical integration is done by providing a 
function that computes the integrand for a given integra- 
tion variable, d, along with the integration limits. The 
function is computed first by solving the equation cf> = 
for y for a given d, then computing the IP terms one at 
a time and summing; each term has a factor which 
corresponds to the product of all the sequences of re- 
flection coefficients mapped by a binary number. A 
numerical integration subroutine is applied which han- 
dles integrands with removable singularities particularly 
well and also provides error estimates. 

5. The current program was tested by considering the 
case when the lower half-space is a vacuum; i.e., the 
case when the structure is a plate. Results obtained are 
identical to the results obtained by the program to com- 
pute the Green's tensor of a plate developed previously 
[12] which had also been checked experimentally. 



12. Computed Results and Discussion 

A FORTRAN program has been developed to numer- 
ically compute the Green's tensor for a layer on a half- 
space with three different bonding conditions between 
the layer and the half-space. The program is written in 
such a way that for given isotropic material parameters, 
maximum observation time, subscript of the component 
of the Green's tensor, number of sampling points, and 
distance between the source and detector, the program 
will compute the displacement at each sampling time. 
Furthermore, the arrivals of various rays are also com- 
puted and identified. The current limitation of the pro- 
gram is that the positions of the source and detector must 
be located on the top surface of the layer. However, the 
program will be modified to be applicable to the case of 
the source and detector located at arbitrary positions in 
the layer. 



Figures 6-27 show the computed results of the com- 
ponents of the Green's tensor and their spatial deriva- 
tives for a plexiglass layer on a glass substrate. They 
were carried out on an IBM compatible personal com- 
puter. The abscissa coordinate, time, is normalized by 
the time required for a shear wave to vertically travel the 
thickness of the layer. The solid curves are the results for 
a welded interface condition, while the dashed curves 
are for a liquid coupled interface condition. The sub- 
scripts ij (11, 12, 13, ..etc.) or ijk (111. ..etc.) indicate the 
response in a specified coordinate direction, i , of a point 
detector located on the top surface of the layer to a point 
force with step function time dependence exerted on the 
top surface in a specified coordinate direction, j . The 
index k denotes the specific spatial derivative direction. 
Physically, the spatial derivative function Gij, k can be 
considered as the displacement in the ('-direction due to 
a differentiated force which is equivalent to a couple or 
a dipole. The distance between the source and detector 
is three times the thickness of the layer (except for Figs. 
7 and 8) and is chosen in such a way that the very large 
Rayleigh wave arrivals are avoided within the given ob- 
servation time, and therefore the details of the early time 
arrivals can be examined. The other reason for choosing 
this distance is to provide a basis of comparison for our 
experiments which we have recently conducted [15] 
with a geometry that is convenient to arrange and carry 
out. It is seen from these figures that the differences in 
the results with different interface conditions are very 
significant. The responses for a delta function time de- 
pendence can be obtained by numerical differentiation 
of the responses mentioned above. The results for differ- 
ent distances are shown in Figs. 6-8. If we compare the 
early time arrivals in Fig. 6 (x = 3) and Fig. 7 (x = 6), 
they are very different only because of the change of the 
distance between the source and detector. This occurs 
because when the distance increases, more head waves 
arrive, and different arrivals change their order of occur- 
rence. The Rayleigh wave arrivals with welded and liq- 
uid coupled interface conditions can be clearly seen 
from Fig. 8 and they are coincident. The early time 
arrivals and their differences between the two different 
interface conditions are totally masked by the amplitude 
of the Rayleigh wave arrivals. Little information about 
the interface can be derived. 

Determination of the Green's function for a structural 
configuration results in a method for determining the 
response of the structure to a temporally varying and 
spatially distributed input loading. The motion of a point 
on the structure can be computed due to, for example, a 
transient contact input force. The transient waveform 
input is the force that would be produced by a point 
source transducer. For an input waveform of interest, the 
motional response is calculated by a point-by-point 
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Fig. 6. Green's function G33. X= 3. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 
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Fig. 7. Green's function G33. X=6. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 
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Fig. 8. Green's function G 33 . X = 1.5. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 
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Fig. 9. Green's function G 2 i- X = 3. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 
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Fig. 10. Green's function G u . X = 3. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 
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Fig. 11. Green's function Gu. X=3. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 
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Fig. 12. Green's function G31. X- 3. (A) The dashed curve is the result for a liquid coupled 
interface. (B) The solid curve is the result for a welded interface. 



0.075 



0.050 



E 

§ 0.025 




-0.025 



-0.050 



P Pp*P p s .p_ 



1.5 



1-" 

pp 



G 331 X=3 -° 



_Pp*S' 



2.0 



Pp*S pp P* pp Sp*S-'" s 



I j 

PS PPPP 



2.5 



3.0 



Normalized Time 

Fig. 13. The spatial derivative of the Green's function G331. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 14. The spatial derivative of the Green's function Gai. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 15. The spatial derivative of the Green's function ftu. X=3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 16. The spatial derivative of the Green's function G131. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 17. The spatial derivative of the Green's function G311. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 



466 



Volume 107, Number 5, September-October 2002 

Journal of Research of the National Institute of Standards and Technology 



G x = 3.0 



o 

z 



0.25 



0.20 



0.15 



0.10 



0.05 



-0.05 



















.'.'' 


.*** 

,...--' 






'^ 


*• T 


,--" 




.,'' 




p' Pp*P-— "" p s *P 

_u t 


Pp*S Ps*E 

L \ 


PPp*PPSp*S s - 


t 

Ipp 


T 

;PS PPPP 



1.5 



2.0 



2.5 



3.0 



Normalized Time 

Fig. 18. The spatial derivative of the Green's function G332- X = 3. (A) The dashed curve is 
the result for a liquid coupled interface. (B) The solid curve is the result for a welded 
interface. 
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Fig. 19. The spatial derivative of the Green's function Gm.. X = 3. (A) The dashed curve is 
the result for a liquid coupled interface. (B) The solid curve is the result for a welded 
interface. 
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Fig. 20. The spatial derivative of the Green's function Gm- X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 21. The spatial derivative of the Green's function G132. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 22. The spatial derivative of the Green's function G312. X=3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 23. The spatial derivative of the Green's function G333. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 24. The spatial derivative of the Green's function G 2 n- X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 25. The spatial derivative of the Green's function Gui-X=3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 26. The spatial derivative of the Green's function G133. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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Fig. 27. The spatial derivative of the Green's function G313. X = 3. (A) The dashed curve is the 
result for a liquid coupled interface. (B) The solid curve is the result for a welded interface. 
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convolution of the time differential of the source wave- 
form with the Green's function component that repre- 
sents the step function response in the motional direction 
of interest. 

The response of a structure due to a damped sinusoid 
input can be easily determined and may be of interest in 
analyzing the response to pulses used for probing mate- 
rials and structures. The response of a layered structure 
interrogated by a pulsed laser beam can also be deter- 
mined by convolving an input waveform with the deriva- 
tive of the components of the Green's tensor that repre- 
sents a dipole source. By comparing the response in the 
case of a vacuum lower half-space with the results of a 
welded solid half-space, the behavior of a large debond 
can be calculated. The large differences in amplitude 
and polarity of the first few head waves and regular 
reflected rays are strong indicators that there is no bond. 
These results and others will be described in a future 
paper. Many other responses of interest for particular 
applications, such as nondestructive evaluation, are en- 
visioned. 



13. Appendix A. The Willis Inversion 
Method 

This appendix directly quotes from Ref. [3]. 
Consider the inversion of transforms of the form 

cc-0< 
-oo-Oj 

exp[-i(£,* B - tat + A^f + K 2 g)] (A.l) 

where F(<u,£) is a homogeneous function of degree —2. 
Each term of the integrand, which represents a general- 
ized ray in the text, is analytic for Im{u>) < and the 
integration of each term vanishes for t < 0. So fit, x), 
i.e., the integration of the sum of the integrands should 
vanish for t < 0. The functions f M 3 and £"3 are homoge- 
neous of degree 1 in a>, £1, and £ 2 ; Ai£ M 3 and A 1^3 have 
imaginary parts. 
Setting 

fl=o»/lfl, i) a = &/\&, (A.2) 

where l£l = (£f + £) m , gives 

oo-oi 00 

/(*,*) = g^lim j dnjis^F(av)jm 

-»-0> lijUl 

exp[-il£l(-/2f + r] a x a + Ai£f + A 2 &" - is)] (A3) 
where £" = i;™(fl, 17). The "convergence factor" e~ 18 is 



inserted to facilitate evaluation of the integrals in Eq. 
(A. 3) by simple quadrature; generalized functions are 
continuous with respect to limiting operations of the 
type introduced, and the limit e — > can be taken at any 
convenient stage. Evaluating the integral with respect to 
l£l gives 



f(f,X): 



cc-0< 



iiji -00-oj 



an 



F(At?) 



(A.4) 



-Ot+ T) a x a + Ai£ M + \i£g - iO/ 1 

Setting 

<j)=-nt+7] a x a + \i^' + \ 2 ^ - is (A.5) 

the integral with respect to Cl is now evaluated by clos- 
ing the contour in the lower half-plane because the inte- 
grand has behavior as 0((T 2 ) for large fl and using 
Cauchy's theorem. This gives 

f(t,x)==^zite— -r^hrf 5r . (a.6) 

477 i-f-! ~ l + Al «.« + A 2fc3,ri , M 

where Q is now taken as the root in the lower half-plane 
of the equation 

-ilt+ J) a x a + Ai^fCAr?) + Kitg{n,T)) = Oi, (A.7) 

if there is such a root. If there is more than one root, then 
the contribution from each root should be included; if 
none, then that term of the integrand is replaced by zero. 
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